//
//---------------------  1D  ---------------------
//
//
// Upwind is not defined for 1D
//

//
//---------------------  2D  ---------------------
//
//       Staggered Mesh for u-vel and v-vel
//
//   0       1       2       3       4       5   
//
//5      >       >       >       >       > 
//       |       |       |       |       |              
//   ^---+---^---+---^---+---^---+---^---+---^  4       Mesh for scalar fields
//       |       |       |       |       | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  3           +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  2           +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//   ^---+---^---+---^---+---^---+---^---+---^  1           0 1   2   3   4 5 
//       |       |       |       |       |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |       |       |       |       |                 x boundary node
//   ^---+---^---+---^---+---^---+---^---+---^  0          > u velocity 
//       |       |       |       |       |                 ^ v velocity
//0      >       >       >       >       >    
//       0       1       2       3       4             
//
//                                                  
//          Volumes for v-velocity
//
//       0       1       2       3       4       
//
//       >       >       >       >       >      5
//       :       :       :       :       :              
//4  ^...+---^---+---^---+---^---+---^---+...^          Mesh for scalar fields
//       |       |       |       |       | 
//       >---o--->---o--->---o--->---o--->      4        5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//3  ^...+...^...+...^...+...^...+...^...+...^              +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//       >---o--->---o--->---o--->---o--->      3           +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//2  ^...+...^...+...^...+...^...+...^...+...^              +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//       >---o--->---o--->---o--->---o--->      2        0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//1  ^...+...^...+...^...+...^...+...^...+...^              0 1   2   3   4 5 
//       |       |       |       |       |                 
//       >---o--->---o--->---o--->---o--->      1          o central node
//       |       |       |       |       |                 x boundary node
//0  ^...+---^---+---^---+---^---+---^---+...^             > u velocity 
//       :       :       :       :       :                 ^ v velocity
//       >       >       >       >       >      0
//
//   0       1       2       3       4       5             
//
//                       
//                  |           |           |           |
//                -->-----o----->-----o----->-----o----->-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |  (i,j+1)  |     :     |
//                  |     ^     |    v_N    |     ^     |   
//                  |     :     |     :     |     :     |
//                  |     :     |     :     |     :     |
//                -->-----o---- 3 -- v_n -- 4 ----o----->--  4 = u(i  , j+1)
//                  |     :     |     :     |     :     |    3 = u(i-1, j+1)
//                  |     :     |     :     |     :     |    2 = u(i  , j  )
//                  |    v_W   u_w   v_P   u_e   v_E    |    1 = u(i-1, j  )
//                  |  (i-1,j)  |   (i,j)   |  (i+1,j)  |
//                  |     :     |     :     |     :     |
//                -->-----o---- 1 -- v_s -- 2 ----o----->-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |     :     |     :     |
//                  |     ^     |    v_S    |     ^     |   
//                  |     :     |  (i,j-1)  |     :     |
//                  |     :     |     :     |     :     |
//                -->-----o----->-----o----->-----o----->--
//                  |           |           |           | 
//                  
//               1           3                       2         4 
//   u_w = ( u(i-1,j) + u(i-1,j+1) ) / 2   u_e = ( u(i,j) + u(i,j+1) ) / 2
//   v_n = ( v(i,j) + v(i,j+1) ) / 2     v_s = ( v(i,j) + v(i,j-1) ) / 2
//               
//---------------------  3D  ---------------------
//
template<class T_number, int Dim>
inline
bool UpwindY<T_number, Dim>::calcCoefficients(const ScalarField &nut) { 
    T_number dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    T_number dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    T_number ce, cw, cn, cs, cf, cb;
    T_number nutinter;
    T_number dxyz_dt = dx * dy * dz / dt;
    T_number RaGaVol = Rayleigh * Gamma * 0.5 * dx * dy * dz;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = ( u(i,j,k) + u(i,j+1,k) ) * 0.5 * dyz;
		cw = ( u(i-1,j,k) + u(i-1,j+1,k) ) * 0.5 * dyz;
		if ( ce > 0 ) ce = 0.0; 
		else          ce = -ce;
		if ( cw <= 0 ) cw = 0.0; 

		cn = ( v(i,j,k) + v(i,j+1,k) ) * 0.5 * dxz;
		cs = ( v(i,j,k) + v(i,j-1,k) ) * 0.5 * dxz;
		if ( cn > 0 ) cn = 0.0;
		else          cn = -cn;
		if ( cs <= 0 ) cs = 0.0; 

		cf = ( w(i,j,k) + w(i,j,k+1) ) * 0.5 * dxy;
		cb = ( w(i-1,j,k) + w(i-1,j,k+1) ) * 0.5 * dxy;
		if ( cf > 0 ) cf = 0.0;
		else          cf = -cf;
		if ( cb <= 0 ) cb = 0.0; 
	    
//
// nut is calculated on center of volumes, therefore, nut
// must be staggered in y direction:	    
		nutinter = 0.5 * ( nut(i,j,k) + nut(i,j+1,k) );

		aE (i,j,k) = (Gamma + nutinter) * dyz_dx + ce;
		aW (i,j,k) = (Gamma + nutinter) * dyz_dx + cw;
		aN (i,j,k) = 2 * (Gamma + nutinter) * dxz_dy + cn;
		aS (i,j,k) = 2 * (Gamma + nutinter) * dxz_dy + cs;
		aF (i,j,k) = (Gamma + nutinter) * dxy_dz + cf;
		aB (i,j,k) = (Gamma + nutinter) * dxy_dz + cb;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) +
		             aN (i,j,k) + aS (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) +
		             dxyz_dt;	    
//		aP (i,j,k) /= alpha;  // under-relaxation
//		+ (ce - cw)  + (cn - cs) + (cf - cb);	    
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.

		sp (i,j,k) = v(i,j,k) * dxyz_dt - 
		    ( p(i,j+1,k) - p(i,j,k) ) * dxz +
		    RaGaVol * ( T(i,j,k) + T(i,j+1,k) ) +
		    nutinter * ( (u(i,j+1,k) - u(i,j,k) - 
				  u(i-1,j+1,k) + u(i-1,j,k)) * dz +
				 (w(i,j+1,k) - w(i,j,k) - 
				  w(i,j+1,k-1) + w(i,j,k-1)) * dx ) +
		    v(i,j,k) * (1-alpha) * aP(i,j,k)/alpha; // under-relaxation
	}    
//    cout << "\n Dim :" << dimension << " function : a_EW_3D "; 

    calc_dv_3D();
    applyBoundaryConditions();
    
    return 1;
}















